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Abstract 

We consider a Frohlich-type Hamiltonian on a hexagonal lattice. Aiming to 
describe nanotubes, we choose this 2-dimensional lattice to be periodic and to 
have a large extension in one {x) direction and a small extension in the other 
(y) direction. We study the existence of solitons in this model using both 
analytical and numerical methods. We find exact solutions of our equations 
and discuss some of their properties. 
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I. INTRODUCTION 



Nanotubes have attracted a large amount of interest ever since they were first discovered 
in 1991 [1]. They can be thought of as carbon cyhnders with a hexagonal grid and are thus 
fuUerene related structures. Their mechanical, thermal, optical and electrical properties 
have been studied in some detail [2]. It was found that most properties depend crucially 
on the diameter, chirality and lengths of the tube. A distortion of the lattice thus affects 
the energy band gap. This distortion of the lattice can be achieved in two different ways: 
a) through an external force like e.g. bending, stretching or twisting [3] or b) through an 
internal excitation, which interacts with the lattice. It is well known that the interaction of 
an excitation such as an amide I- vibration in biopolymcrs or an electron (in the case of the 
Frohlich Hamiltonian) with a lattice whose distortion is initially caused by the excitation 
results in the creation of a localised state which, in what follows, we refer to as a soliton. 
Such a soliton was first introduced by Davydov [4] in the 1970s to explain the dispersion 
free energy transport in biopolymers (see also [5] for further details). 

Recently, a Frohlich Hamitonian was studied on a two-dimensional, discrete, quadratic 
lattice [6-8]. In [6,7], the existence of localised states was studied numerically and it was 
found that their properties depend crucially on the electron-phonon coupling constant. An 
analytical study confirmed these results [8] - by showing that in the continuum limit the set 
of discrete equations reduces to a modified non-linear Schrodinger (MNLS) equation which 
has an additional term resulting from the discreteness of the lattice. Although a soliton of 
the basic NLS is unstable the extra term was shown to stabilise it for appropriate choices of 
the coupling constant. 

In this paper, we extend the study of Refs. [6-8] to the case of a hexagonal, periodic lattice 
with a large extension in x and a small extension in y directions. We study the resultant 
equations both analytically and numerically. In Section II we present the Hamiltonian and 
the equations of motion. In Section III we discuss various properties of the equations in 
the stationary limit and we demonstrate the existence of an exact solution of the discrete 
equations. In this limit, we can thus replace the full system of equations by a modified 
discrete non-linear Schrodinger (DNLS) equation. In Section IV, we present and compare 
our numerical results for the continuous MNLS equation, for the full system of equations 
and for the modified DNLS equation. Most of our numerical simulations are performed for 
the (5, 5) armchair tube. Such nanotube was discussed in a recent paper by Liu ct al. in 
the context of nanorings [9]. We also discuss briefiy how nanotubes with different diameter 
and chirality could be constructed in our model. 

II. THE HAMILTONIAN AND EQUATIONS OF MOTION 

A. Hamiltonian 

The Hamiltonian H of our model is a sum of four sums which result from the special 
features of the hexagonal grid. V^j^ denotes the electron field on the i-th, j-th lattice side, 
while Uij and Vij are the displacements of the i-th, j-th lattice point from equilibrium in 
the X direction and y directions, respectively: 
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jx is the electron field self-interaction coupling, Cx couples the electron field to the displace- 
ment fields u and v and kx is the self-coupling of the displacement fields. 



B. Equations of motion 

We can easily derive the equations of motion from our Hamiltonian H. As an example 
we give the equations for i = 1 + ik. The discrete Schrodinger equation for the ipij field 
thus becomes: 



dib ■ ■ 

^{E + w)i^i^j - 2jx + A-i,j + A+i,j-i) 



while the equations for the displacement fields and arc given by: 



+ ^ (siV'i-ijf - IV'i+ij+iP - IV^i+ij-iT) 



and 



df 



= kx (3vjj - Vj+ij+i - Vi^i^j - -yj+ij-i) 
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The equations then read: 
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III. STATIONARY LIMIT 

In the stationary limit, we have: 

+ xjjij[Ui+ij+i + C/,+i,,_i - 2C/,_i,,- + V3(l^,+i,,-+i - y^+ij-i)] = (11) 
with A = £'o + Wo — 6 and 

3Uij - Ui+ij+i - Ui^ij - C/i+ij-i = -| (2|V'i-ijr - \fpi+i,j+i\^ - iV'i+ij-iH (12) 
3Vij - Vi+ij+i - Vi^ij - Vi+ij^i = ^ (iVi+ij+iT - lA+ij-iH (13) 

with ^ = 

A. Discrete equation 

In contrast to the square grid, we find that the discrete equations of the hexagonal grid 
in the stationary limit do have an exact solution. We can thus replace the system of coupled 
equations (3)- (5) by just one modified DNLS equation. We again look at the case i — l + 4k 
for which we have : 

mUij = I (21^.-1,/ - iV'mj+il' - iV'w-iH 
where A(l)C/jj = C/j+ij+i + Ui-ij + C/j+ij-i — SUij. Analogously, we have: 

A(1)V^..- = ^ (l^mj-ir - IV'mj+ir) 

where A{l)Vij = V^+ij+i + ^-ij + V^+ij-i — SV^j. Next we note that for the three nearest 
neighbours we have similar relations, namely: 
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A(i)c/,+i,,_i = I (iV',,,-2r + iV'M-r - 2iVi+2j-ir) , 
A(i)t/.-i,, = I (iV'.-2,,-ir + iV'.-2,,+ir - 

for the [/ field and 

A(l)ym,.>i = ^ (lV'.,f - lV'.,.+2r) , 
A(l)T^m.-i = ^(lV'..-2r-|V'd') 

for the y field. 

Defining: Z„ = t/j+ij+i + t^i+ij-i - 2t/i_ij + \/3{Vi+ij+i - l^+ij-i) (i.e. the lattice 
terms in (11)) we find that the following discrete equation holds 

A(1)Z,= 

9 (eiV'ijf - \ipi,j+2\^ - \'ipi+2,j+i\'^ - i^M-2r - iV'i+2j-ir - \A-2,j-i\'^ - iV'i-2j+ir) 

The right hand side of (14) is a 7-point Laplacian A(2)|'0jjp, thus we find: 

A(l)Z„ = -^A(2)|V^,,r (14) 
It is easy to see that one possible solution of this equation is of the form: 

Za = -g (iV'i+ij+il^ + + iV'i-ijf + SiV'ijf ) . (15) 

This is quite remarkable since on a square lattice a similar equation has no simple solution. 
Inserting (15) into (11) we have: 

Xijjij + 2 {3ipij - iJi+ij+i - ^i-ij - ^i+ij-i) 
- + + l^i-ijf + 3|^ijf ] = 

or 

Xjjjij - 2A{l)^ij - gij,,, (A(l)|^,jf + 6|^,,f ) = 0. (16) 
This equation constitutes our discrete nonlinear Schrodinger (DNLS) equation. 
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B. Continuum limit 



Next we look at the continuum limit of (16). To do this we introduce the following 
expansions: 

dijj dil) 1 O'^ip d^^l^ 1 O'^ijj 



2 



dx dy 2 dx"^ dxdy 2 dy 



and 



. - . - . ... (18) 

where for i = 1 + 4/c and i = 3 + 4A'; we have = 1/2, 5a;_ = 1, while for i = 2 + Ak and 
i = 4 + 4A; we have = 1, = 1/2. Moreover, 8y — -\/3/2. Inserting this into (16), we 
obtain: 

Xi, - ^Aij - -gij Qa|^|2 + 6|^|2^ = (19) 

or, equivalently, 

~XiP + Ail; + Agip(^\iP\^ + ^A\iP\^^^0 (20) 

We thus have, in analogy to what was found in [6-8], a modified nonlinear Schrodinger 
(MNLS) equation with an extra term, which can stabilise the soliton: 

M + AtI> + Ag^ fl^p + lA\tl;A = 0. (21) 



Following [6-8] we see that the conserved energy in this case is given by: 

E = I (^1 - 2^|V'|' + f (VIV'I^)^) dxdy . (22) 

Approximating the soliton solution of (21) by a Gaussian of the form ijj(x,y) — 
exp(— ^(x^ + y"^)) we find that the value of ac that minimises the energy is 

- 2(1 - |) (23) 
and thus we have an estimate of the critical g, namely, (jcr tt. 
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IV. NUMERICAL RESULTS 



A. Continuous, modified non-Unear Schrodinger (MNLS) equation 

First, we have considered the continuous, modified non-hnear Schrodinger equation (21). 
For this, we have taken the radially symmetric ansatz: 

V^(r,t) = e-*i?(r), (24) 

and put this into (21). Then we have solved the ordinary difi^erential equation (ODE) using a 
collocation method for the boundary value ODEs [10] and choosing the boundary conditions: 

^^|,=o = , i?(r = oo) = 0. (25) 
or 

Our results are shown in Fig. 1. For a fixed g we have determined the value of the function 
R{r) at the origin, R{r = 0), as well as the value of a for which the norm of the solution 
is equal to unity. As can be seen from Fig. 1, we have found that for a critical value of 
9 — 9cr ^ 2.94 the value of -R(O) tends to zero. Since our construction is such that the 
maximum of the solution is located at r = 0, the value of the height of the solution tends 
to zero and thus the solution ceases to exist. We thus find that the critical value of g from 
our numerical construction agrees with the upper-bound value obtained from the variational 
approach based on the Gaussian, namely, gcr — tt. 

Note that the solutions cease to exist when the value of a tends to zero. Since wc can 
interpret a as the frequency of an internal rotation, the solutions apparently cease to exist 
when there is no internal rotation. This can be compared to the so-called "Q-balls" which are 
non-topological solitons characterized also by a complex scalar field [11] of the form similar 
to (24). For them its is known [11] that there exist upper and lower positive- valued bounds 
on the frequency of the internal rotation in order for Q-balls to be stable. In comparison, 
our solutions exist for all values of a > 0. This is probably due to the fact that while the 
dynamical part of our action is similar to that for Q-balls, we have a extra term involving 
derivatives as compared to an "ordinary" (/)*'-potcntial in the case of Q-balls. 



B. Discrete equations 

1. Full system of equations 

For our numerical study of the full equations (3)-(5) we have found it convenient to 
"squeeze" the lattice as indicated in Figs. 2a und 2b. The Hamiltonian and the corresponding 
equations are given in the Appendix. For our numerical calculations, we have used mainly 
a periodic grid with A^i = 160 and N2 = 20. Wc have in addition chosen the boundary 
conditions such that the fields at {i = 0,j) are identified with those at {i = imax,j)- Thus 
the type of nanotube we are studying here is a (5,5) armchair tube which is metallic [2]. 
Nanotubes can also be semiconducting and we make a brief comment about the possibility 
of constructing semiconducting tubes in our model in the last paragraph of this section. 

In this work, the C-C bond length, 0.1415 nm [12], is normalised to unity. Therefore, 
the tube diameter is 0.6756 nm. Tubes with different diameter can also be constructed in 
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our model. We discuss this together with different chirahties in the last paragraph of this 
section. 

As starting configuration we have used an exponential-like excitation ip^j extended typ- 
ically over the lattice points i = 78 — 83 and j — 3 — 7 with the lattice at equilibrium 
everywhere, i.e. Uij — and Vij — for all We are mainly interested in the existence 
of solitons and their dependence on the value of the coupling constant c^. Wc have set 
jx = kx = 1, M = 20 and E = 0.142312. The main goal of this work is to study the 
dependence on c^. So the exact values of and play a minor role. Hence, we have set 
them to one. The choice of M = 20 is a reflection of the physical fact that the mass of the 
carbon atom is approximately 20 • 10~'^^g. 

To absorb the energy thus allowing the initial configuration to evolve into the stationary 
solutions of (3)- (5), i.e. of (11)-(13) we have additionally introduced damping terms i'^^^ 
and i^-^, respectively, into the equations (4) and (5). We have typically chosen u = 
0.25 — 0.75. For this choice of the coupling constants, we have performed several numerical 
calculations using a 4th order Runge-Kutta method for simulating the time evolution. We 
have found that sohtons exist in this system for Cx ^ 20. For larger values of Cx, the 
soliton forms very quickly, while decreasing Cx the time increases at which a soliton forms. 
This is of course due to the weaker coupling between the dynamics of the lattice itself 
and the excitation. For Cx = 19, we have waited until t ~ 8000 and haven't found a 
soliton. Moreover, in all cases we have found only little displacement of the lattice from 
the equihbrium. We have found that at the location of the soliton the lattice becomes 
squeezed (i.e. the lattice sites move towards the sites at which the soliton is located). This 
is demonstrated in Fig. 3 for Cx = 25, where we show the lattice distortion after t = 4000. 
The point at which the centre of the soliton is located doesn't move, while the sites in its 
close neighbourhood all move towards the centre of the soliton. 

We have also studied the effects of perturbations of the solitons. We have found that 
after perturbing the soliton we obtain a new solution with a different height of the soliton 
maximum. Even after introducing a perturbation which keeps the maximal height fixed, 
the new solution differs from the starting one. We thus come to the conclusion that the 
full system of equations (3)-(5) has a large number of solutions for each choice of couphng 
constants. We believe that a conserved quantity exists in this system which picks out the 
specific solution. However, so far we have not been able to determine this conserved quantity. 



2. Modified, discrete non-linear Schrddinger (DNLS) equation 

In addition to the full system of equations, we have also studied the dynamical analog 
of equation (16). Using a similar starting configuration with ip^^j being exponential and 
nonzero over i = 78 — 83 and j = 3 — 7, we have determined the value of g for which a 
soliton exists. Our results are shown in Fig. 4 where we present the height of the soliton's 
maximum {ipip*)max as function of g. We find that the value of g at which the soliton 
disappears gcr ~ 2.295. The height of the soliton at this critical coupling is {ipi,jipl^j)max ~ 
0.227. Our numerical study of the continuous MNLS equation gave us gcr ~ 2.94, while the 
analytic study led to gcr — t^- Both values are not a bad approximation for the value found 
numerically for the discrete equation. 
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To test the independence of our results from the form of the initial settings, we have 
used a different starting configuration with two exponential- like excitations being located at 
i = 78 — 83, j = 3 — 7 and i = 138 — 143 and j = 13 — 17, respectively. We have found that for 
values ^ ~ 3, the results agree. For both types of initial configurations, the minimal energy 
configuration corresponds to one soliton. However, having said this, the time to reach this 
minimal energy configurations is significantly smaller for the initial configuration with one 
excitation than for that with two (typically one order of magnitude smaller). We have also 
tested our results as to the dependence on the size of the grid. For this, we have chosen 
two excitations on three different grid sizes: a) a grid with A^i = 160, N2 — 20 and two 
exponential excitations extended over i = 78 — 83, j = 3 — 7 and i — 138 — 143, j — 13 — 17, 
respectively, b) a grid with A^i = 320 and = 40 with the excitations located at the 
same places than in a), and finally c) a grid with A^i = 60, N2 = 10 and two exponential 
excitations extended over i = 18 — 23, j = 2 — 4 and i = 38 — 43, j = 7 — 9, respectively. We 
have found that for ^ = 3, the results of cases a) and c) agree. For the case a) the sohton 
forms at t ^ 300, while for the case c) it forms a,t t ^ 100. This is not suprising since in the 
case c), the two excitations are located nearer to each other than in the case a). To test the 
dependence on the actual lattice size we have compared the cases a) and b). We have found 
that the larger the lattice the longer it takes for the soliton to form. For ^ = 3 a soliton 
forms after t ^ 300 in the case a), while for b) it forms at t > 700. We have thus found 
that, in comparison with the case of the full system of equations, the solutions of the DNLS 
equation are unique for each choice of the coupling constant. 



3. Comparison of results 

Since we have found that in the stationary limit the full system of equations can be 
replaced by a DNLS equation, the minimal energy solutions we have obtained for both 
types of equations should be in agreement. 

Comparing the two systems, we see that the value g is given in terms of the coupling 
constants of the full system by: 

2 

^^QMjX ^^^^ 
which, for the choice of coupling constants we have used in our numerical simulations, gives: 

9 = — ■ (27) 
^90 ^ ^ 

Thus a critical value of ~ 20 would imply (jcr ~ 4.4. First, we remark that the values of 
the critical electron-phonon coupling we obtained from all our simulations (including those 
for the continuous MNLS equation) are of the same order of magnitude. However, there is 
a slight discrepancy between the results for the full system and the DNLS equation. We 
believe that this is due to the fact that there might exist additional terms A in (15) for 
which A(l)^ = and/or A(2)A — 0. These terms would then appear in (16) and would 
change the comparison of the solutions. However, it is difficult to determine these additional 
terms and so this is left as a future work [13]. 
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4- Tubes with different diameter, chirality and lengths 

Since most of our results are for a (5,5) armchair tube and since it is well known that the 
properties of nanotubes depend strongly on the diameter, chirality and lengths of the tube, 
we will discuss briefly how different tubes could be constructed in our model. We haven't 
constructed these tubes yet, but we aim to do so in a future pubhcation in which we intend 
to extend our approach to a more reahstic three dimensional model [13]. 

Labelling the first carbon atom in y direction by j = 0, we have chosen j„iax such that 
it is divisable by 4. Thus, the length of the tube in the y-direction is ly = \imax- Since we 
identify the fields labelled by (i = 0, j) with those at (i = imax,j), the diameter of the tube 
is given by d = ly/n. Thus increasing/decreasing jmax by 4n, n = 1, 2, 3..., we can construct 
armchair nanotubes with diameters d = -^{jmax ±4n). Similarly, we can construct longer 
tubes by increasing the number of atoms in the x direction. 

As far as chirality is concerned, there are two things to modify in our model in order to 
be able to construct tubes with different chirality. One is to change the number of points 
in the y direction so that jmax is non-divisable by 4. The other is to adjust the periodic 
boundary conditions in the y-dircction appropriately. If we e.g. choose jmax = 18, we have 
to identify the fields at {i = 0,j) with those at {i = imax,j + !)• This then would give us a 
(5,4) nanotube which would be semi-conducting. 

In this work, we have concentrated our attention on the existence of localised structures. 
These structures extend over large parts of our grids, but arc negligible at boundaries. 
Hence, we expect these to hold for systems with different boundary conditions, i.e. different 
chiralities. 

V. CONCLUSIONS 

Motivated by a large amount of research done in the area of nanotubes, we have studied 
solitons on a 2-dimensional hexagonal lattice. We have chosen our lattice to be periodic in 
both the X and y directions and to be of large extention in one (the x) direction and of small 
extention in the other (the y) direction. In the stationary limit, we have found that the full 
system of equations in which the electron excitation is coupled to the displacement fields of 
the lattice can be replaced by a modified discrete non-linear Schrodinger (DNLS) equation. 
This discovery of an exact solution of the full system of equations is remarkable since for 
the similar quadratic lattice such a simple solution does not exist. 

In our numerical studies we have mainly concentrated our attention on determining the 
value of the critical phonon-electron coupling constant. For the DNLS we have found that 
unique solutions exist and that the value of the critical coupling is in good agreement with 
both the analytically and numerically values found for the continuous analog of the DNLS. 
For the full system of equations, we believe that a large number of solutions exist for each 
choice of the coupling constants and that a conserved quantity exists in the system. The 
critical value of the electron-phonon coupling is of the same order of magnitude as in the 
case of the DNLS; however, we believe that this small discrepancy results from the fact 
that possible "boundary" terms appear when replacing the full system by the DNLS. These 
boundary terms are terms which are annihilated by either the 4-point Laplacian A(l) and/or 
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by the 7-point Laplacian A(2). To find these terms is non-trivial and since this seems an 
interesting topic by itself, we leave this as a future work [13]. 

Finally let us mention that a possible extension of the results given here would involve 
the study of the corresponding three dimensional equations and/or of the influence of exter- 
nal forces. 

Acknowledgements BH was supported by an EPSRC grant. 
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VI. APPENDIX 



A. Hamiltonian and Equations of motion for the numerical studies 



To simplify the numerical construction of the solutions we have squeezed the lattice as 
indicated in Fig.s 2a, 2b. This reduces the memory requirements and so speeds up the 
calculations. The Hamiltonian H"^ for the numerical construction thus takes the form 



N2 



JVi_ 



(c c 



^1 1 



j=l 1^=0 



N2 -T- 



.=1 i=i 



with the phonon energy : 



+ 9^ E E (^^[(^y - ""i-ij+i)^ + (% - f^i-ij+i)^]) 



(28) 
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1 ■ ^ ^ 



The equations of motion are then given by 
• For z = 1 + 4fc, A; = 1,2...: 

dipi ■ 



dt 



- 3^ (2V'i-ijV'r-ij - A+iji^i+ij - V'i+ij-iV'r+ij-i) 



—kx (Svij — Vi+ij — Vi-ij — Vi+ij-i) 



For i = 2 + 4A;, A; = 1,2... 



ih- 



dt 



{E + W)'il)i^j - 2jx {i/^i+ij + + V'i-ij+i) 



d'^u 



dt^ 



d^v. 



df^ 



+ 



(V'i-ij+iV'tij+i - V'i-ijV'tij) 



For i = 3 + 4A;, A; = 1,2...: 

Oijji ■ 



dt 
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(37) 



— i'^Vij — Vi+lj — Vi-ij — Vj+ij+i) 



(38) 



For i ^ 4: + 4:k,k = 1,2...: 



dt 



j{ui^i,j + -Ui-ij-i - 2iii+ij) - - (39) 



(40) 



+ 



(41) 
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FIG. 1. The results obtained with the ansatz i/j{r, 0) = e*"*i?(r) to solve the continuous modified 
NLS equation. In this figure, the value of the function R{r) at the origin, i?(0), which in our 
construction is equal to the maximum of R, Rmaxi is shown as function of g. The values of a are 
also shown. The solutions shown have norm one. 




FIG. 2a. The hexagonal lattice is shown. The arrows indicate the method of "squeezing" the 
lattice for the numerical evaluation. 
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FIG. 2b. The "squeezed" hexagonal lattice for the numerical construction. The Hamiltonian 
ffn corresponding to this lattice is given in the Appendix of the paper. 




45 46 47 48 49 50 51 52 53 54 55 

X 

FIG. 3. The distortion of the lattice close to the location of the soliton is shown. The squares 
indicate the undistorted lattice, while the circles indicate the distorted lattice after t = 4000; 
Cx = 25. The corresponding soliton's maximum {tpip*)max ~ 0.6145. 
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FIG. 4. The height of the sohton's maximum {'4'i,j'4^ij)max is shown as a function of the pa- 
rameter g. 
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